/*
 * 蒙特卡洛积分：求单位半圆面积
 */

#include<iostream>
using namespace std;
#include<cmath>

double f(double x){
	return sin(x);
}

void mc_integ(double (*func)(double x), double a, double b, double & NIf, double & error, int n){
	
	int i;
	double x, f, sum_f = 0, sum_f2 = 0, ave_f, ave_f2, I1, I2;

	for(int i=0;i<n;i++){
		x = (double)rand()/RAND_MAX*(b-a) + a;//取随机数 x 在 [-1,1] 均匀分布
		f = func(x);
		sum_f += f;
		sum_f2 += f*f;
	}
	ave_f = sum_f / n;
	ave_f2 = sum_f2 / n;

	NIf = (b-a)*ave_f;
	error = 8 * sqrt( (ave_f2 - ave_f * ave_f)/n );
	// 4 sigma 置信度为 99.993666%
}

int main(){

	double NIf, error;
	int n=1E4;
	mc_integ(f, 0, M_PI, NIf, error, n);
	cout<<" 4 sigma accuracy 99.993666%: If = "<< NIf <<" +- "<<error<<endl;
	return 0;
}
